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We investigate growth on vicinal surfaces by molecular beam epitaxy mak- 
ing use of a generalized Burton-Cabrera-Frank model. Our primary aim is 
to propose and implement a novel analytical program based on a perturba- 
tive solution of the non-linear equations describing the coupled adatom and 
dimer kinetics. These equations are considered as originating from a fully mi- 
\Q \ croscopic description that allows the step boundary conditions to be directly 

' formulated in terms of the sticking coefficients at each step. As an example, 

■ we study the importance of diffusion barriers for adatoms hopping down de- 

. scending steps (Schwoebel effect) during growth and post-growth equilibration 
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O . I. INTRODUCTION 

O 

_ Growth on vicinal surfaces, with a misorientation angle of few degrees, is of profound 

^ ' importance as a method of preparing high-quality epitaxial layers and heterostructures.lll 
■ This is also frequently used to obtain insight into the fundamental kinetics of molecular beam 
epitaxy (MBE) growth.ii The renewed interest in the analytical model of Burton, Cabrera, 
and Franki (BCF) describing epitaxial growth on vicinal surfaces is thus not surprising. 
In the BCF model, a crystal grows by step propagation due to incorporation of atoms 
deposited onto the terraces. No interactions between the atoms are considered, and it is 
also supposed that adatom concentrations at the step edges are at equihbrium and the 
effect of the step movement can be neglected. This model has been recently generalized for 
applications pertaining to growth by MBE.i^ The far-from-equilibrium character of MBE 
has led to lifting; all of the above restrictions and taking into accoimt the influence of the 
step movement ,0 deviations from local equilibrium at the step edges,El and, finally and most 
importantly, effects related to lateral interactions of adatoms on the terraces, in particular 
formation of adatom islands.0^ 

Our primary concern in this paper is to propose and carry out an analytical program for 
solving the generalized BCF equations including non-linear terms due to adatom interac- 
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tions. The non-linear structure of these equations has been the major obstacle to obtaining 
analytical solutions and for this reason only numerical solutionsIrEl have been reported so 
far. The steady state solution of these equations^ was previously considered through the 
use of a microscopic formulation based on adatom and dimer position-velocity distribution 
functions. f{x, v), J-'{x, v), instead of the macroscopic densities m{x), M{x) more commonly 
usedjU'ErEi and we will implicitly use this formulation here as well. However, as discussed 
below, here this distinction is only operationally significant in the treatment of the boundary 
conditions which we are able to describe in terms of the step sticking coefficients. 

We will assume that the growth conditions (substrate temperature, growth rate, etc.) are 
such that step flow dominates so that only small islands (dimers) have to be considered, but 
the step movement is still sufficiently slow so that the steps can be considered stationary in 
solving the boundary value problem.Q'l These restrictions are imposed only as simplifications 
and their removal would increase the complexity and decrease the clarity of our results. All 
other assumptions introduced in this paper will be clearly noted and assessed in the following 
sections. While it is our belief that some of these could be eliminated, our main goal here is 
to obtain a closed form analytical solution that is based on sound physical and mathematical 
reasoning. 

The outline of this paper is as follows. The necessary equations and notation are pre- 
sented in Sec. |11 A| where the distinction between the microscopic and macroscopic formu- 
lations is made explicit. Our procedure for the solution of these equations then follows in 
Sec. |II B| and this is implemented in Sec. |T|. As a demonstration of our results, we study 
the effects of barriers to hopping down descending step edges (Schwoebel effect) in Sec. 
Sec. |V] is a brief summary of the results of the paper. 



II. THEORETICAL CONSIDERATIONS 
A. Basic Equations 

We consider a stepped crystal substrate on which adatoms are uniformly deposited at a 
rate Fa^ where F is the flux (m~^s~^) and a the lattice constant (see Fig. |I|). As indicated 
in the introductory section, the temperature and the flux are supposed to be such that only 
dimer formation has to be considered but the steps are still sufficiently slow and can be taken 
as stationary in solving the boundary value problem.0 Equations describing the evolution of 
the adatom and dimer densities have been presented and used by several authors. These 
same equations can also be obtained from the microscopic description used by us earlier in 
treating the steady state proble m0 if we also assume that the dimers are immobile and that 
the adatom current is given by Pick's law. The last assumption, which is an integral part of 
the macroscopic description, is equivalent to neglecting an initial stage in which the current 
changes over the time scale tq (see below) from zero at the time t = to the Pick's law 
value. This should not have a significant qualitative or quantitative impact on our results 
and therefore seems justified at this juncture of the theoretical development. 

Denoting the adatom and dimer densities by m{x,t), M{x,t) to indicate that thev are 
the result from averaging the distribution functions /, over velocity, we then haveuQ'Hy 
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^ = + F + 27M - (8/r)m - 2Tm' (1) 

^ = (4/r)m + rm2-7M (2) 
< a; < L , t > , 

where D is the diffusion coefficient, r = 7 = Tq"^ exp(— i^^/Zc^T) with the charac- 

teristic time for adatom surface diffusion I/tq = 2D /o? = pq exp(— E'd/Zcb?'), z/q being the 
attempt frequency, En the dimer binding energy and Eo the surface diffusion barrier for a 
free adatom, and T is the temperature and ks is the Boltzmann's constant. Also, we have 
taken the number of dimer-forming sites around a free adatom as four.illi 

The terms on the right-hand side of Eq. (P represent (from left to right): diffusion of 
free adatoms, deposition of an atom on a surface site with all nearest-neighbor positions 
unoccupied, dimer break-up, dimer formation due to direct collision of an atom from the 
beam with a surface adatom, and dimer formation when two diffusing adatoms meet each 
other (cf. Fig. 0). The prefactor F = aa? /t^ with a the capture efficiency which we sub- 
sequently set to unity,0 combines a geometrical factor cra^ and a relaxation time that we 
have denoted as x^. In previous numerical studie^Hi Tj. = tq has been used; here we will use 
Tr = (FD)^^/"^ which is a considerably larger quantity. It can be shownEl that Tr = (FD)^^^'^ 
is in quantitative agreement with Eqs. (0) while r,. = tq is not. Separate justifications 
based on kinetic arguments have also been given in the literature.0 

The initial conditions for Eqs. ([^), are m(x, 0) = M{x,0) = 0, i.e. the beam is 
initiated at t = at which time the terraces are unoccupied. So far there is no difference 
between our description. Ens. (|T|), (0) with the above initial conditions, and the macroscopic 
descriptions cited earlier .llHl However, in specifying the boundary conditions the advantages 
of the microscopic theory become clear and a major difference between the two approaches 
emerges. At the level of approximation in the microscopic theoryll^llll leading to the above 
kinetic equations, the distribution function / describing the adatoms is given as 

f{x,v,t) = |mi(x,t)e(t;) + m2(a;,t)[l-0(t;)]|(27ra)-^/2g-«V2a (3) 

where G(f) is the standard Heaviside function and a = kBT (the adatom mass is expressed 
in units in which it is unity). This leads to 

m= {mi + m2)/2 , j = = (a/27r)^/^(mi - 1712) (4) 

It follows from Eqs. (0) and (0) that mi is the density of adatoms moving with positive 
velocity on the terrace and 1712 the density of those moving with negative velocity. This 
level of detail is particular to the microscopic description and allows us to prescribe natural 
boundary conditions. The most general of these, which we use here, is given in terms of the 
step sticking coefficient Si, i = 0,L for x = 0,L, respectively (Fig. 0). For asymmetric step 
conditions at x = 0,L 

So = K(0) - mi(0)] /m2(0) , Sl = [m,{L) - m2(L)] /m,{L) , (5) 
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with < 5*1 < 1. The finite parameters 5*^ describe the full range of step conditions from 
blocking {Si = 0) to total absorption {Si = 1). In addition to the conceptual advantage 
in using the microscopic boundary condition there is also an operational advantage. From 
Eq. (D we find 

din 

— him = — — , x = L (6a) 

ox 

dm 

hom = — — , x = (6b) 
ox 

with hi = [Si/D{2 - Si)] (2a/7r)i/2. Thus, for the full range of step kinetics the boundary 
condition is of the Neumann (or radiation) type and a single solution depending on the St 
describes all possible sticking conditions. Eqs. (H), and (||) together with the initial 
conditions define the problem to be considered in mathematical terms. Note that there 
is no boundary condition for M required for the case considered of slow steps and fixed 
dimers. In the steady state an exact solution for these equations was obtainedli^ despite 
their non-linearity; here we cannot expect to be so fortunate and therefore we have to look 
for approximate means of solution. 



B. Plan of Solution 

We propose to solve the problem defined above by tailoring a procedure of solution 
to our intuitive expectation of how the physical process evolves. Initially the terraces are 
unoccupied, adatoms are deposited and begin to diffuse to the steps and, as their density 
increases, dimers begin to be formed and, finally, the latter may begin to disassociate. 
Thus, three separate evolutionary stages, most likely not sharply defined, are suggested. 
One mathematical interpretation of this is to consider the terms in Eq. (|1|) representing 
dimer effects to be of higher order in an approximation scheme; M also will be of higher 
order. Writing m = m^^^+m^^^+ . . . , M = M^^^ + M^'^^ + . . . and re- writing Eq. (|ip we have 

^ = ifi(m,Af) (2') 

{R{m, M) represents all the terms on the right-hand side of Eq. (^ which are not given 
explicitly) so that our procedure for solution is implemented by considering the following 
equations, 

^..^_.(„.o),0, (S) 
^ = ifl(m(»),MW), ... (9) 



with, e.g. 
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(10) 



The connection between the preceding physical and mathematical descriptions becomes 
clearer if we operate on Eq. (|l]) with d/dt and use Eq. (^) to obtain 



d_ 

di 



dm ^ d'^m ^ 
dt dx"^ 



-1 



dm ^ d'^m ^ 
dt dx"^ 
dR2{m) 



d_ 
di 



8m/T + 2Tm' 



-Ri{m) - 



dt 



(11) 
(I'O 



The first term on the right-hand side is the contribution due to dimer break-up, indi- 
cated by the 7 prefactor, while the second is the contribution due to dimer formation. Since 
i?^(m(o))=0 by virtue of Eq. (0) and 9i?2(m^°^)/5M0 we see that the mathematical descrip- 
tion given by Eqs. (j^), (||) reflects our conjectured physical picture in which dimer break-up 
effects are of higher order than dimer formation effects. In lowest order the right-hand side 
of Eq. ([TT|), which contains only dimer effects, is set equal to zero and the left-hand side used 
to determine m^^\ which identically recovers Eq. since dF /dt = 0, dm^^^x , 0) / dt = F , 
d^m (°)(x,0)/9x2 = so that 



d_ 

dt 



dt 



-D- 



dx'^ 
dm^^\x, t) 
dt 



d_ 
dt 



dt 



dx^ 



, and 



^ d'm('\x,t) ^ ^ 
dx"^ 



(12a) 
(12b) 



In what follows we will directly use Eqs. (|^ and (§) to determine m^°\ m^^\ 

The preceding discussion contains a rationale, based on physical reasoning, together 
with a supporting mathematical argument, for implementing a perturbative approximation 
of Eqs. (1'), (2'). Our reasoning here is admittedly heuristic, however given the intractable 
nature of reaction-diffusion equations in general (a bounded domain is the most difficult 
case), we believe that the approach outlined above is a reasonable first step in obtaining a 
useful analytical solution for the generalized BCF equations in the regime near step flow 
where dimer effects begin to play a role. It has been shown! that in this regime the moving 
boundary effect can still be neglected so that this omission here is not an issue. It is 
possible that the procedure outlined above could be formalized through the use of a scaling 
argument leading to the introduction of one or more parameters of smallness; we have not 
attempted to do this here since the underlying physical rationale is so intuitively appealing 
and transparent. In the next section we obtain explicit results for m^^\ m^^\ and M(i) by 
implementing the program described above. 



III. SOLUTIONS 

A. Lowest Order; m^*^) 

As discussed above, we write m = m'^'^'>+m'^^'>+ . . . and determine m^^\ m^^^ from Eqs. (^, 
together with the boundary conditions specified by Eq. (|TUp. These equations can be 
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solved by standard methods;EJ we first consider the equation for m^^^ , Eq. (|^) for the general 
case where the beam F is turned off at t = t*. The exact solution for t<t* is 



r Li 

:,t) = J2 (F/Dbl) (l - e-""^*) x) / x')dx' 

n=l ^ 

where the eigenf unctions K{bn,x) areEl 

K{bn,x) = V2[bnCosbnX + hosinbnx] \{bl + hl){L + hLibl + hl)-^ + ho) 



-1/2 



(13) 



(14) 



and the eigenvalues bn are given by the positive roots of b{hQ + Hl) = (6^ — /io/^-l) tanfoL. 
The terms without the exponential factor can be directly summed but the tedious algebraic 
manipulation required can be avoided by noting that these terms are the non-decaying 
part of the solution, i.e. the steady state solution of Eq. (|^. In order to simplify the 
remaining sum and obtain a result that is both easier to use and also makes the next level 
of approximation tractable we will replace 6„, in the exponential term only, by a lumped 
eigenvalue, 6.0 In Sec. [Ill D| we indicate a procedure for determining b in terms of the system 
parameters through the imposition of a consistency constraint. 
Making the above simplification we then have 



m 



(0) 



[x, t) = {F/2DAo) (-Aox^ + h^B^x + B^) To{t) 



(15) 



where 



Bo 



ho + hL + hohiL 
2hoL + hohLL'^ 
1 - e-^'^* , t<t* 



t > t* 



The above results can be non-dimensionalized (see Appendix 
convenient to retain the dimensional form shown above. 



(16a) 
(16b) 
(16c) 
(16d) 

but for now it is most 



B. First Order: Adatom Density m^^-* 
The equation for m^^'' is given from Eq. (^) as 

^ - = -(8/r)m(o) - 2r(mW)2 = R,{m^'^) (17) 

This equation is a more general version of Eq. (|^), i.e. a diffusion equation with a source 
term that here is not constant but depends on both space and time. The formal solution is 



CO ,.f 

m« = ye-'"^*ir(6„,a;) / dt'e'"^*' / dx' K{bn,x')R2{m^''\x' ,t')) (18) 
r^i ^0 -^0 



We denote the contribution from that part of linear in m*^'') as m*^^^^ so that m^^^ = 
_l_j7^{i2). fu];ther we note that m'^^^^ will be identically zero when the beam is turned 
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off. Using Eq. (|T5|) for m^^^ and again replacing 6„ by b in the exponential terms after the 
t' integration is carried out, we find 



AF 
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X* + 



h^B 
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X + 



where 



(19) 



All 
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'Aoh 






12 


^11 (t) 







+ h^Bo - BoL 



(20a) 
(20b) 

(20c) 



The contribution from that part of R2 non-linear in m^^\ denoted previously as m^^'^\ is 
obtained from Eq. (|18]). Making use again of the simplification introduced earlier we obtain 



m 



(12), 



x,t) = -2T 



2DA 
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■X — X H X + 
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(21) 
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bhLj 12 



ihoB,2y 



2A12-B12 — 



6Ai2Bi2ho 



h, 



{hoBu? 



'2A12B12 



hoBf2 

hr 



+ 



^12 



~2b''Dte-'"''' 



t < t 



-b'^D(2t-t*) 



t > f 



(22a) 

(22b) 
(22c) 

(22d) 



C. First Order: Dimer Density M^^^ 



In first approximation the dimer density is found from Eq. (|^); although this equation 
includes disassociation effects these do not contribute directly to the solution which is con- 
sistent with the logic underlying the approximation procedure we are using.lli The initial 
condition for Eq. (^ is M^^\x,0) = M{x,0) = 0, so that the direct contribution of the 
disassociation term 'jM^^^ vanishes and we obtain 



2 Jo 



[X, t ) 



(23) 
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so that M(^) is given directly from Eq. ([TsD following a simple integration. If we write 
^(0) _ fn(^^ (^x)J-'Q{t) then, in notation identical to that used for indicating the separate 
contributions to m*^^-' 



(24) 



^2i(t)=7-'(l-e-^*)-(7-&'^) 



-l/-b'^Dt 



--yt\ 



t < t 



t > f 



(25a) 
(25b) 



M(i2)=r(mW(a;))V22(t) 



(26) 



^22(t) = 7-'(l - e-^*) - 2(7 - 62l^)-i(e 



l/~2b'^Dt 



--yt\ 



t < t 



{-f-2b'D)-\e 
^-i(e-7(*-**) _ e-7t) _ _ 62^)-i(e-^'^*^e-^(*-**) - e"^*) 
(7 - 2b^Dy^e'^' + (7 - 2b^Dy^ 



t > r 



(27a) 



(27b) 



D. Determination of the Lumped Eigenvalue b 

The major approximation we have made is replacing 6„ by b in Eq. (|T3p which is the 
basis for Eqs. (|15D , (|19|), and (^). It would be possible to use a second lumped eigen- 



value in the calculations used in obtaining the latter two results, but this would appear to 
over-complicate matters at this point. What is important for our present purposes is the 
prescription of a rule that fixes b in terms of the system parameters. The basis for such a 
rule is a consistency requirement that the value of the coverage as determined from Eqs. (P, 
(I) and d^) be identical with that found directly from their solution; we show how this can 
be done using the zeroth-order approximation. 

If we determine the coverage by directly integrating the solution, m^^\ we find 

e(t) = Tn^^\x,t)dx = 
Jo 

{FiyUD) UL{ho + Hl) + h^hLL^ + 12) {h^ + Hl + h^hLiy^ Mt) (28) 



Using Eqs. (|l|) and (^, after integrating the former we find 

dQ 



-/iLm(°) (L, t) - /iom(°) (0, t)\ (29) 
from which it follows that 

e(t) = {FL/b^D)J^o{t) (30) 

and then 

= (L/12) (4L(/io + Hl) + hohLL^ + 12) {ho + + hohLL)'^ (31) 
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E. Discussion 



Before considering some numerical consequences of the above solutions, which we do in 
the next section, it is best to add here a few final comments pertaining to the procedure 
followed in obtaining these. The approximation procedure is ad hoc and follows from physical 
reasoning supported by the consistency of the resulting mathematical formulation. Although 
we have only considered the first approximation beyond the step-flow solution, in principle 
higher-order corrections could be obtained. It appears to be more useful to consider instead 
additional effects due to moving steps and larger islands in the present context. The former 
of these appears to be tractable, and we hope to be able to report on this in a future 
study. Finally, we want to emphasize that the procedure used here is distinct from standard 
methods based on linearization or a formal small parameter expansion and in our opinion 
captures the essence of physical processes that are taking place better than either of these 
methods. 



IV. EXAMPLE APPLICATION: SCHWOEBEL EFFECT DURING GROWTH 

AND POST GROWTH RECOVERY 

A. Surface roughness and its evolution during growth and post— growth equilibration 

of the surface 

The roughness of growing surfaces has become one of the main topics of interest in the 
study of surface processes. One reason for this is the evident importance of minimizing 
surface roughness for applications, e.g. for the deposition of narrow layered structures 
(quantum wells, lateral superlattices, magnetic multilayers, etc.) required for the fabrication 
of novel devices such as quantum well lasers or modulation-doped field effect transistors. 
The performance of such devices depends critically on the morphological quality of epitaxial 
layers and interfaces. 

Another reason for the intense interest in surface roughness is purely theoretical. The 
roughness of growing surfaces has been observed to exhibit asymptotic dynamical scaling be- 
havior which has led to the classification of growth models into various universality classes. 
In a majority of the theoretical studies, the focus has been on the evolution of the sur- 
face morphology due to fluctuations in the incoming flux of particles and surface diffusion. 
However, recently the influence of additional activation barriers to adatom hopping between 
layers (the Schwoebel effect (SE), Ref. [T^J20) on the developing surface roughness has been 



studied by Villain0 and identified as a cause of an intrinsic growth instability on nomi- 
nally flat surfaces. On vicinal surfaces, the additional step-edge barriers to hopping down 
descending steps stabilize growth and lead to equalization of terrace widths.SS 

From a practical point of view, the growth interruption technique@ has become a useful 
method for obtaining atomically abrupt interfaces during growth by MBE. In this technique, 
growth is briefly interrupted by stopping the cation (but not the anion) flux. After the 
growth is stopped, the surface smoothes which causes the recovery of the specular-beam 
intensity of the commonly used monitoring probe, reflection high-energy electron-diffraction 



9 



(RHEED), to its initial pre-growth value. It has been found experimentally that the recovery 
of the intensity I{t) can be fitted by an empirical expressionc3 

/(t)=Ao + Aie-*/"i+A2e-*/"2 ^ t>t* (32) 

where Aq, Ai, and A2 are constants, and ri and T2 are the time constants for the fast and 
slow stages of recovery, respectively. The time constant of the initial stage of recovery, ri, 
has been shown to have a strong systematic dependence both on the substrate temperature 
and on the point in layer completion at which growth is interrupted.!! 

The recovery process is also of profound theoretical importance. The equilibration of 
a surface after a period of growth is a more discriminating test of theoretical models than 
growth alone. During growth under typical conditions, the maximum time scale is set 
by the deposition flux and processes occurring over longer time scales are "frozen out". 
However, during recovery, these processes can come into play even though their effect during 
growth can be safely omitted.@ Indeed, Monte Carlo (MC) simulation studies which yielded 
very good results when directly compared to observed behavior during growth, appeared 
to be unable to reproduce the recovery curves, in particular at low temperatures, where no 
temperature dependence of ti was observed in simulations.! Similarly, the detailed study 
of the recovery process by Vvedensky and ClarkeB using the same MC model showed that 
in their simulations there was no relevant dependence of ri on the point where growth was 
interrupted. 

Peng and Whelanll suggested that inter-layer transport is of key importance for the 
correct reproduction of the behavior of ti. However, large values of ri observed in experiment 
can be explained only if there is some mechanism hindering inter-laver transport. The 
obvious choice for such a mechanism is the SE, and recent simulationsH^ have shown that it 
indeed results in slowing down the recovery process. However, it has been also found that the 
SE alone cannot explain the observed behavior and an additional process must be introduced 
in which a freshly arrived atom is incorporated at a site with the highest number of nearest 
neighbors in a certain region around the initial site of incidence.!!! This incorporation scheme 
results in smoother edges of both pre-existing steps and adatom islands. An independent 
proof exists that such a change in the model leads to a better reproduction of the post- 
deposition recovery. In Ref. the authors had to consider detachment of atoms from steps 
and break-up of islands composed of up to 10 atoms in order to obtain good qualitative 
agreement with the experimentally observed behavior. Significantly, they assumed a shape 
for each island that is as close as possible to a square {i.e., very compact and with smooth 
edges). Similarly, the break-up of the most stable configurations has been found to be a 
crucial step for the post-growth relaxation of the surface in Ref. 

MC simulations which included both the SE and the smoothing due to the incorpo- 
ration have led to quantitative agreement with the experimental dataS which is a strong 
circumstantial evidence for the existence of the Schwoebel barriers on GaAs(OOl) surface. 
Additional evidence has been provided in a recent paper by Johnson et a/.il who observed 
(using scanning tunneling and atomic force microscopies) growth instability and formation 
of large mounds on a nominally fiat GaAs(OOl) surface. Finally, a recent analysisil of the 
RHEED data of Ref. |^ demonstrated that sticking coefficients at ascending and descending 
steps on GaAs(OOl) surfaces are indeed strongly asymmetric providing further support for 
the existence of the step-edge barriers on this surface. 
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In light of the above, it would be interesting to study the influence of the SE using the 
analytical results presented above and taking advantage of the natural treatment of the 
boundary conditions. We will use the parameters used for simulations of GaAs(OOl) growth 
and try to make a comparison to some of the experimental results of Ref. 



B. Adatom and dimer concentration profiles and densities, the surface step density, 

and their time evolution 

In this section, we discuss adatom and dimer concentration proflles, and the time evo- 
lution of the total densities of adatoms and dimers and of the surface step density, during 
deposition and after its cessation. We study both the case with and without the SE {i.e. with 
different and equal sticking coefficients at boundaries of the terraces, respectively). Note 
that the non-dimensionalized form of the equations has been used in numerical evaluation 
of the equations, cf. Sec. |ITT| and Appendix 0. 

In our calculations, we proceeded with "deposition" until the saturation of the total 
density of adatoms and dimers on the terraces (i.e. a steady-state situation) has been 
achieved, see Fig. ^, and then we interrupted it. The terrace width L has been chosen to 
be 20 lattice constants. The sticking coefficient at the ascending step {Sq in Fig. |ip has 
been estimated based on the dimer binding energy E^- (assuming straight steps similarly to 
Myers-Beaghton and Vvedensky in Ref. ||, Sec. IVB), as Si = l—exp{—EN/kBT), whereas the 
sticking coefficient at the descending step {Sl in Fig. 0) was estimated using the additional 
step-edge barrier height Eb as well. Si = exp{—EB / ksT) —exp^—E^ / ksT) . The values of 
parameters E£, = 1.54 eV (which determines the diffusion coefficient D, see the paragraph 
following Eq. (^ in Sec. |II A] ), E^ = 0.231 eV, and Eb = 0.175 eV were taken from the 



results of the fitting procedure in Ref. |2^. The attempt frequency z/q has been calculated 
from i>o = 2kBT/h, where h is the Planck's constant. The value of the flux F was 0.20 ML/s, 
and the values of the substrate temperature T are given in the flgures. 

Adatom and dimer concentration proflles shown in Fig. ^ have been calculated at the 
time t* when the deposition was stopped (in the saturation region). The upper two panels 
(a and b) of Fig. |^ are for the case without the SE, whereas the lower two (c and d) for the 
case with the SE. The panels a and c show both the total adatom density m (solid lines) and 
the adatom density in the lowest order of approximation, m^'^^ (dotted lines), cf. Eqs. (|15|) 
and ([H). 

The proflles for the case with the SE are asymmetric and exhibit a maximum near the 
reflecting boundary. Notice that signiflcantly higher adatom and dimer densities build up on 
the terraces due to the presence of the Schwoebel barriers. It is also interesting to mention 
that whereas the densities with and without the SE are of approximately the same order of 
magnitude just after the growth is stopped (Fig. |^) they differ by many orders of magnitude 
after the equilibration due to a larger time constant for this process in the case with the SE 
(see below). Note also that concentration proflles do not vanish at boundaries, x = 0, 1 (in 
contrast to those in Ref. |7|,|30D due to the microscopic treatment of the boundary conditions 
in our approach, cf. Sec. |II A . 



In Fig. 1^, the time evolution of the total densities of adatoms and dimers on terraces is 
shown. The strong temperature dependence of the densities and of the difference between m 
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and m^^^ can be again observed. However, the most significant difference is a rather drastic 
change of the time for the decay of the adatom and dimer densities due to the SE. This 
confirms that the inclusion of the barriers to hopping at step edges does lead to a slower 
relaxation of the surface after growth is terminated, a conclusion reached previously in the 
MC study mentioned above.0 

Finally, in Fig. ^, the evolution of the surface step density calculated in a pseudo-two- 
dimensional approximation (similar to that used in Ref. |^) is shown. For the calculation 
of this quantity, we supposed that every free adatom and dimer contributes 4 and 6 steps, 
respectively. The panel a shows the results without the SE, whereas the panel c those 
with the SE. The experimental results at similar temperatures taken from Ref. § (RHEED 
specular-beam intensity during the post-growth recovery on a GaAs(OOl) vicinal surface 
misoriented by 2° toward [010] direction) are shown in the panel b. Notice the order of 
magnitude difference in the time scale which is discussed in the following section. 



C. Discussion 



The above numerical results are based on evaluation of the preceding analytical results 
and agree with earlier numerical studies. The advantage of our approach (besides the 
important fact we were able to obtain analytic formulae) is the microscopic level of the de- 
scription of the growth process and resulting natural treatment of boundary conditions. We 
have demonstrated than even a relatively small Schwoebel barrier has a profound influence 
on behavior during growth and, in particular, during the post-growth recovery. Another 
detail worth mentioning is that had we used tq instead for the time constant in the pref- 
actor r in Eq. (|^) (see the discussion in Sec. [11 A] ), the contribution m^^^^ to the adatom 
density would be underestimated by several orders of magnitude. Let us now discuss several 
problems related to our results. 

As expected, the difference between the adatom concentration profiles obtained using 
formulae for m and m^^^ (Fig. |^ a, c) is not important under conditions near the step-flow 
growth mode (the biggest difference being observed near the maxima), but increases as 
the temperature is decreased (cf. similar results in Fig. 1 of Ref. ^). However, this does 
not mean that the dimer (and larger islands) formation is not important. It was shown by 
Stoyanovill and Myers-Beaghton and Vvedensky0 that adatom interactions lower the adatom 
diffusivity by several orders of magnitude and, if not taken into account properly (as it is the 
case when the Einstein relation for surface diffusion = 2Dt is use J), lead to an estimate 
of the critical temperature Tc (at which growth starts to proceed via step flow) which is 
incorrect by up to two hundred kelvin.lll In other words, whereas the concentration profiles 
can be calculated in the lowest order of approximation near T^, the critical temperature 
itself must be determined taking higher-order effects into account. 

The time constants for the post-growth recovery of the adatom and dimer densities are an 
order of magnitude smaller than observed in experiment (cf. Fig. |^ and particularly Fig. 
This is not surprising given the fact that in our model, the effect of break-up of larger 
islands is not taken into account. It was convincingly demonstrated by Myers-Beaghton 
and Vvedenskyi and Kenny et a/.0 how significant these processes are for the correct repro- 
duction of the recovery curves, and this process was also shown to be important in the MC 
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study mentioned abovecj where it was found that step and island edges smoothed by the in- 
corporation mechanism of the incoming atoms (see above) act in concert with the step-edge 
barriers to slow down the recovery process. Also the fact that one-dimensional diffusion in 
our model leads to more rapid changes of the surface adatom population as compared to the 
two-dimensional case adds to the difference in the speed of surface smoothing. 

Finally, in the framework of our simple model, we cannot observe oscillations of the den- 
sity of adatoms, dimers, or the surface step density (or oscillations of related quantities like 
the specular-beam intensity of some diffraction probe calculated in the kinematic approxi- 
mation). Such oscillations are not experimentally observed near the step-flow growth mode 
where our model is, strictly speaking, only supposed to be valid. However, it is possible 
to obtain oscillations at lower temperatures (or higher fluxes) even in this class of models 
provided the step movement is taken into account. □ As we already mentioned, the inclusion 
of this effect into our perturbative scheme seems to be possible (cf. Sec. |I|). 

V. SUMMARY 

In summary, we have investigated epitaxial growth on vicinal (stepped) surfaces using 
an analytical program based on a perturbative solution of the generalized Burton-Cabrera- 
Frank equations. The non-standard scheme followed by us is motivated by physical insight 
into the growth process. This approach allowed us to obtain closed-form analytical solutions 
for the adatom and dimer densities while using a microscopic description of the boundary 
conditions. The scheme is open to further improvements and inclusion of other effects such 
as the movement of the pre-existing steps. 

We applied the analytical results to the study of the influence of the Schwoebel effect 
(additional barriers to inter-layer hopping at step edges) on the evolution of adatom and 
dimer populations and of the surface step density during growth and post-growth equi- 
libration of the surface and compared them (where it was possible) to the experimental 
results. The Schwoebel effect was shown to have a significant impact on the behavior during 
growth and, in particular, during the post-growth recovery, but it has also become apparent 
that larger islands and their decay have to be taken into account if the recovery curves 
obtained experimentally and in simulations are to be reproduced. In order to obtain a more 
realistic description of the growth process at lower temperatures (including oscillations of 
step density or kinematic intensity), additional processes must be considered such as the 
above-mentioned step movement. 
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APPENDIX A: DIMENSIONLESS VARIABLES 



It is convenient for the purposes of assessing the effects of various parameters to represent 
the results obtained in dimensionless form. Introducing x' — xL and {hi)' = ahi and using 
the relationships D — a'^/2To, a/To — y/2a which provide reasonable estimates, we then find, 
e.g. for m^°^ 

a2m(°) = {To/rAo)[-Ao{L/afx^ + hoBo{L/a)x + Bo]To{t) (Al) 

where we have dropped the primes for notational convenience. Now < x < 1 and 

= [S,/{2 - Si)] (2/0F) (A2) 

Ao = ho + hi + hQhi{L/a) (A3) 

Bo = hi{L/a)^ + 2{L/a) (A4) 

and 

l/Db^ = (L/a) (ro/6A) [A{L/a){ho + h) + hoh{L/a)^ + 12] (A5) 

Similar results follow for m^^^ and M^^^ but because they are so lengthy we do not include 
them here. 
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FIG. 1. Schematic picture of the processes considered in our model. Atoms are deposited (1) at 
a rate F onto the substrate with a lattice constant a held at a temperature T. They migrate (2) 
with the diffusion coefficient D and attach either to the pre-existing steps (with probabilities given 
by the sticking coefficients Sq and Sl) or form dimers (3), (4), which can subsequently disassociate 
(5). Two ways of a dimer creation, by a deposition of an incoming atom onto a nearest-neighbor 
site of a surface adatom (3) or by an encounter of two migrating adatoms (4), are shown. 



FIG. 2. Adatom (a, c) and dimer (b, d) concentration profiles immediately after growth is 
stopped at the steady-state region. The panels a, b and c, d show the results obtained without 
and with the Schwoebel effect, respectively. See text for more details. 



FIG. 3. Time evolution of adatom (a, c) and dimer (b, d) densities. The panels a, b and c, d 
show the results obtained without and with the Schwoebel effect, respectively. 



FIG. 4. Time evolution of the surface step density (panels a, c) and the RHEED specular beam 
intensity (b). The panels a and c show the results obtained without and with the Schwoebel effect, 
respectively. The direction of vertical axis in the plots of the step density is reversed to allow 
for comparison with the RHEED specular-beam intensity (see Refs. |,26), i.e. the step density 
increases downwards. 
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